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RUNGE-KUTTA METHODS 

FOR LINEAR ORDINARY DIFFERENTIAL EQUATIONS 


D.W. ZENGG AND T.T. CHISHOLM 


Abstract 

Three new Runge-Kutta methods are presented for numerical integration of systems of linear 
inhomogeneous ordinary differential equations (ODEs) with constant coefficients. Such ODEs 
arise in the numerical solution of the partial differential equations governing linear wave 
phenomena. The restriction to linear ODEs with constant coefficients reduces the number of 
conditions which the coefficients of the Runge-Kutta method must satisfy. This freedom is used 
to develop methods which are more efficient than conventional Runge-Kutta methods. A 
fourth-order method is presented which uses only two memory locations per dependent variable, 
while the classical fourth-order Runge-Kutta method uses three. This method is an excellent 
choice for simulations of linear wave phenomena if memory is a primary concern. In addition, 
fifth- and sixth-order methods are presented which require five and six stages, respectively, one 
fewer than their conventional counterparts, and are therefore more efficient. These methods are 
an excellent option for use with high-order spatial discretizations. 


Introduction 

We consider the numerical integration of large linear inhomogenous systems of ordinary 
differential equations in the form 

= Au-g(/) 0) 

dt 

where A is an M by M matrix whose elements depend on neither u nor t, and u and g(t) are 
vectors of length M. Such essentially autonomous systems arise in the numerical solution of 
partial differential equations (PDEs) governing linear wave phenomena after application of a 
spatial discretization such as a finite-difference, finite-volume, or finite-element method. 
Examples of such PDEs are the linearized Euler equations governing acoustic waves and the 
Maxwell equations governing electromagnetic waves. The elements of A depend on the PDE 
and the spatial discretization. The inhomogeneous term g(f) is associated with either a source 
term or the boundary conditions. In the context of wave propagation, the system of ODEs is 



often mildly stiff with the eigenvalues of A typically lying near the imaginary axis. 

The system of ODEs arising from the application of a spatial discretization to a system of 
PDEs can be very large, especially in three-dimensional simulations. Consequently, the 
constraints on the methods used for integrating these systems are somewhat different from those 
which have driven much of the development of numerical methods for initial value problems. 
Due to their high accuracy and modest memory requirements, explicit Runge-Kutta methods 
have become popular for simulations of wave phenomena [5,6,7,15,17]. Third- and fourth-order 
methods requiring only two memory locations per dependent variable are particularly useful 
[3,13,14]. This property is easily achieved by a third-order Runge-Kutta method [14], but an 
additional stage is required for a fourth-order method [3]. Since the primary cost of the 
integration is in the evaluation of the derivative function, and each stage requires a function 
evaluation, the additional stage represents a significant increase in expense. For the same reason, 
error checking is generally not performed when solving very large systems of ODEs arising from 
the discretization of PDEs. 

There have been several attempts to develop efficient methods for integrating linear systems 
of ODEs [4,9,10,1 1]. The basic premise of these methods is that the major cost in evaluating the 
derivative function is in forming the matrix A and the vector g(r). In the application considered 
here, the simulation of linear wave phenomena, the matrix A is never explicitly formed or 
stored. Hence the methods previously proposed for linear systems are not appropriate for this 
application. 

It is well known that a Runge-Kutta method with p stages has an order of accuracy not 
exceeding p [1,2]. For p< 4, methods of order p can be derived with p stages. However, fifth- 
and sixth-order methods require at least six and seven stages, respectively. Nine stages are 
required for seventh-order accuracy and eleven for eighth-order accuracy [1]. Since the cost for 
our application is roughly proportional to the number of stages, this represents a significant 
limitation of higher-order Runge-Kutta methods. 

Several authors have considered various approximations to reduce the number of stages and 
the storage requirements of high-order Runge-Kutta methods. Shanks [12] was able to develop 
schemes with a reduced number of stages by requiring only that the accuracy conditions be 
approximately satisfied. Zingg et al. [16,17] propose methods with low storage requirements 
which are of high order for linear homogeneous ODEs but second-order otherwise. A similar 
idea was proposed previously by Lorenz [8]. 
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In this paper, we develop Runge-Kutta methods specifically for linear ODEs with constant 
coefficients. By removing the constraints imposed by nonlinearity in the derivative function, 
high-order Runge-Kutta methods can be derived which are more efficient in some respect than 
the classical methods. In the next section, we present a fourth-order method which requires less 
memory than the classical fourth-order Runge-Kutta method. \Ve then present fifth- and sixth- 
order methods requiring fewer derivative function evaluations per time step than fifth- and 
sixth-order Runge-Kutta methods applicable to nonlinear problems. 


General Form of an Explicit Runge-Kutta Method 

Without loss of generality, we consider the following scalar ODE: 

du r \ 

at 

A general p-stage explicit Runge-Kutta method can be written as 
*1 = f(t n , u n ) 

i-l 

k, = f (t„+Cih, Un+h^aijkj) i = 2, • • • ,p 

j = i 

p 

Wfl+l = Un + hJ^bjki 
i = l 

where h = At is the time step, t„ = nh , and u„ is an approximation to u ( t n ). 


(2) 


(3) 


Low-Storage Fourth-Order Method 

We consider first the case p = 4. With the constraints 

c 2 = «2I 

c ?> = a 31 + a 32 ( 4 ) 

c 4 = a \\ + <242 + a 43 

there remain ten parameters. For fourth-order accuracy, there are eight conditions which must 
be satisfied. Four of these arise even for linear homogeneous constant-coefficient ODEs. A 
further three conditions must be met if the ODEs are inhomogeneous. The final condition is 
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associated with non-constant coefficients or nonlinearity. Therefore, fourth-order Runge-Kutta 
methods are a two-parameter family of which the classical method is a particular choice. 

If we restrict our attention to linear constant-coefficient ODEs, the number of conditions is 
reduced to seven. These are 

4 

X*. = 1 

J=1 

X c,A = 1 n 

i=2 

c 2 a 32^3 +^4( c 2°42 + c 3 a 43) = ^ 

C 2 a 32^43^4 = 1/24 (5) 

= 1/3 

1=2 

= 1/4 

1=2 

^3^2 a 32 + ^4(C2 a 42 + c 3°43) = 1/12 

The reduction in the number of conditions to be satisfied does not permit us to reduce the 
number of stages. However, we can obtain reduced storage requirements. 

Following the approach of Wray [14], the requirement that only two memory locations be 
used imposes the following three constraints: 

b { = a 4 1 = a 3l (6) 

b 2 = a 42 

With these constraints, only two memory locations are required for both the dependent variable 
and the value of the time derivative. Hence the method requires minimal storage even when 
compact or spectral methods are used for the spatial discretization. With the memory locations 
denoted A and B, the method proceeds as follows. 

1. Initially, u„ is stored in A, and B is empty. 

2. The term k \ = /(/„,«„ ) is evaluated and stored in B. 
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3. The quantity u n + ha 3 [ k j , is calculated and stored in A. 

4. The quantity u n + ha 2 \ k \ is calculated and stored in B. 

5. The term k 2 = /(f n +C2^’ u n +ha 21 k \ ) is evaluated and stored in B. 

6. The contents of the two memory locations are linearly combined to form 

u„ + h (a 3 1 k \ +a 22 k 2 )> which is stored in B. 

7. With a 41 =031, another linear combination gives u n + h (a 41 k \+a^ 2). which is stored 
in A. 

8. The term k 3 =/[f„+c 3 /i, w n +/i(a 3 i*i+«32*2)] is evaluated and stored in B. 

9. The contents of the two memory locations are linearly combined to form 

u n + h (041 k i +fl42^2+ a 43^3)> which is stored in B. 

10 . With^i =041 and b 2 = «42> another linear combination gives u n + h (fci/c 1 +t 2 ^2 + ^3^3)> 
which is stored in A. 

11. The term £4 = fUn+cth, M„+/i(fl4i/ci+042^2 +fl 43^3)] * s evaluated and stored in B. 

12. The contents of the two memory locations are linearly combined to form u n+ \ 

With the additional constraints imposed by the low-storage requirement, we are left with 
seven parameters to satisfy the seven conditions given in eq. (5). Although this system may 
possess more than one solution, the only solution we have found is 

a 2 ] =c 2 = 0.69631521002413, c 3 = 0.29441651741, 

c 4 = 0.82502163765, b x =a 4 , =a 3 i =0.07801567728325, 

a 32 = 0.2 1 6400840 1 3679, b 2 = a 42 = 0.04708870 117112, 

a 43 = 0.6999 1 725920066, b 3 = 0.47982272993855, 

b 4 = 0.39507289160708 


Five-Stage Fifth-Order Method 

For the case p=5, we have, in addition to the constraints given in eq. (4), the following 
condition: 
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c 5 = a 51 + a 52 + °53 + a 54 


(7) 


Consequently, adding the fifth stage has produced five additional parameters for a total of fifteen. 
The coefficients must satisfy the following eleven conditions in order to produce fifth-order 
accuracy for linear constant-coefficient ODEs: 

= 1 
i=l 

5 

XcA = 1/2 
1=2 

c 2 a 32^3 + ^4( c 2 a 42 +c 3 a 43) + ^5( c 2 a 52 +c 3 a 53 +c 4°54) = 1/6 

f2 a 32 a 43^4 + ^5[ c 2( fl 42 a 54 + fl 32 a 53) + c 3 fl 43 a 54l = 1/24 

c 2 fl 32 fl 43 a 54^5 = 1/120 (8) 

IV, ? = 1/3 

i=2 

IV? = 1/4 
1=2 

^3C2 a 32 + ^4( c 2 fl 42 + c 3 a 43) + ^5( c 2 a 52 + c 3 a 53 + c 4 a 54) = 1/12 

IV? = 1/5 

1=2 

(^>44142 + ^3^32 + ^5 a 52) c 2 (^5 a 53 ^4 a 43) c 3 ^5 a 54 c 4 = 1/20 

^5t«54(^42 c 2 +a 43 c 3) + a 53 a 32 c 2] + 2 , 4 a 43 a 32 c 2 = 1/60 

Thus a four-parameter family of solutions is obtained. Several different criteria can be 
applied in order to choose a method from this family. The following values have been found by 
minimizing the L 2 norm of a vector containing the coefficients of the method: 
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a 2 l =c 2 = 0.21, c 3 = 0.43, c 4 =0.68, c 5 =0.85, 

a 32 = 0.474 1 85463659 15, a 42 = 0. 1 3437223603429, 
o 52 = 0.26302355344001, b 2 = 0.16574368303091, 
o 43 =0.57068167533284, o 53 = 0.10434139625551, 

b 3 =0.41041645692809, o 54 = 0.39377303853165, 

b 4 =-0.04092124960122, b 5 = 0.37240141 154501 
b, =0.09235969809721 

with 

0 3 i=C3 - 032, 041 =C 4 -a 42 -«43. 

0 5 1 =C 5 -O 52 -«53 ~ a 54 


Six-Stage Sixth-Order Method 

With p= 6, the following condition must be satisfied in addition to the constraints given in 
eqs. (4) and (7): 

C 6 =fl 61 +a 62 +fl 63 +a 64 +a 65 

Therefore, there remain twenty-one free coefficients. The requirement of sixth-order accuracy 
for linear constant-coefficient ODEs produces the following sixteen conditions: 

6 

= 1 

1=1 

= 1/2 

i=2 

C 2 a 32^3 + b 4 (c 2^42 + C 3 a 43 ) + b 5 (c 2 a 52 + c 3 a 53 + c 4 a 54) 

+ b(,( c 2 a 62 + c 3 a 63 + c 4 a 64 + c 5 a 65) = ^ 

C 2 «32«43^4 + ^5[ c 2( a 42 a 54 + «32 a 53) + ^3^43054] 

+ ^6[«65( fl 54 c 4 +^53^3 +«52C2) + «64(043C3 + a 42C2) + a 63«32^2] = 1/24 
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a 65[ a 54(°43 c 3 + a 42 c 2) + a 53 a 32 c 2] + «64 a 43 a 32 c 2 


+ ^5 a 54 a 43 a 32 c 2 


c 2°32 a 43 a 54 a 65^6 ~ 1/720 


1/120 

( 10 ) 


T,b iC } = 1/3 

i= 2 

EV/ = 1/ 4 

i=2 

b 3 c\ a 32 + ^4( c 2 a 42 + ^3 a 43) + ^s( c 2 a 52 + c 3 a 53 + c 4 a 54) 

+ b fc(a 65C5 + fl 64^4 +«63 c 3 +a 62 c 2) = ^2 

£*>,4 = 1/5 

i=2 

1 3 

(^4 fl 42 + ^3^32 + 65052 + ^6 a 62) c 2 + (^5 a 53 + 64043 + ^6 a 63) c 3 

+ 5 « 54 + ^6 fl 64) c 4 + ^6^65^5 = 1/20 

^6l«65(«52C2 +^53^3 + O 54 C 4 ) + <*64(«42C2 + a 43 C 3 ) + a 63 a 32 cl] 
+ £»5[a54(a42^2 + a 43 c 3 ) + a 53 a 32 c 2] + ^4 a 43 a 32 c 2 = */60 

£6 f c? = 1/6 


^6( a 62 c 2 +a 63 c 3 + a 64 c 4 + fl 65 c 5) 

+ b 5 (a 52 C 2 +a 53 C 3 +a 54 C4) + b 4 (a 4 2C2 +a 43 C3) + b 3 a 32 C2 = 1/30 
c\ [^6( a 65 a 52 + a 64 a 42 + 063 ^ 32 ) + 65(054042 + a 53 a 32 ) + b 4 a 43 a 32 ] 

+ c 3 [bf,(a^ 3 a 33 + a^4a 43 ) + b 3 a 34 a 43 ] + c 4 b^a^ 3 a 34 = 1/120 

^6t a 65 a 54( a 43 c 3 + a 42 c 2) + a 64 a 43 a 32 c 2 ] + ^5 a 54 a 43 a 32 c 2 = 1/360 

Using the same criterion as for the fifth-order method, the following coefficients have been 
chosen from the five-parameter family of solutions to the above conditions (again possibly 
nonunique): 
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with 


021 =c 2 =0.15, c 3 = 0.36, c 4 = 0.57, 

c 5 = 0.75, c 6 = 0.90, 032 =0.45818181818182, 


a 42 =0.09769454545455, 
a 62 =0.20874226393025, 
a 43 =0.48766666666667, 
a 63 =0.12686271445897, 
a 54 =0.44703799502007, 
6 4 =0.35718962665957, 
65 =-0.00478351095633, 
, = 0.03850905269576 


a 52 =0.10861879806510, 
b 2 =0.24971305394585, 
a 53 =0.04655817933320, 
=0.11278150363005, 
a M =0.02734417934727, 
a 65 =0.37591957583530, 
£> 6 = 0.24659027402511 


a 31 — C3 — 032, 04] — c 4 0 42 0 43> 

051 = C 5 - a 52 -053 “0 54. 


0 61 = c 6 _ «62 “ 0 63 “064 _a 65 


Stability Contours 

The stability contours of the three new methods are shown in Fig. 1. Satisfaction of the first 
four conditions in eq. (5) ensures that the new fourth-order method has the same stability contour 
as the classical fourth-order Runge-Kutta method. Similarly, the stability contours of the five- 
stage fifth-order method and the six-stage sixth-order method are uniquely defined and do not 
depend on which members of the respective families are selected. 

Although the stable regions of the fifth- and sixth-order methods are somewhat larger than 
that of the fourth-order method, the increase is not sufficient to compensate for the cost of the 
additional stages. Therefore, the fourth-order method is a better choice if the time step is limited 
by stability considerations. The stable regions of the fifth- and sixth-order methods do not 
include the imaginary axis. Systems with pure imaginary eigenvalues are obtained when central 
differencing is applied to the spatial derivatives in partial differential equations governing wave 



propagation phenomena with no physical dissipation, in the absence of boundary conditions. 
However, Zingg et al. [17] have demonstrated that by adding a small amount of numerical 
dissipation to the spatial discretization, stable schemes can be obtained using such methods. The 
amount of dissipation required is sufficiently low that the overall accuracy of the scheme is not 
compromised. The stability contour of the method successfully used in [7] for simulations of the 
propagation and scattering of electromagnetic waves is identical to that of the present sixth-order 
method. 


Fourier Error Analysis 


Using Fourier analysis we can determine the errors produced by an integration method when 
applied to a linear homogeneous ODE. Since our interest is in wave propagation, we consider a 
scalar ODE of the form 


du 

~dt 


i(SM 


( 11 ) 


where to is a real constant. The Runge-Kutta methods developed here produce a solution in the 
form 


«„ = <fuQ ( 12 ) 

where 

o = £ 77 ( 1 ( 0 )* (13) 

k= 0 * • 

and p is the number of stages. The local amplitude and phase errors are determined from o as 
follows 


er n = I at - 1 


er p = — 


tan (a,/a r ) 

co/i 


+ 1 


where a r and a, denote the real and imaginary parts of a. 


(14) 

(15) 


Figs. 2 and 3 show the local amplitude and phase errors produced by the three new methods. 
In order to account for the number of stages, the errors are plotted versus coh/p. Hence the errors 
shown are for approximately equal computational effort. Since the time step is thus proportional 
to p, the amplitude error shown is |a| 1 /p -l. The figures show that each increase in the order of 
the method produces an increase in accuracy even though the extra work has been accounted for. 



Hence the fifth- and sixth-order methods can be more efficient than the fourth-order method if a 
sufficiently accurate spatial discretization is used. 


An Example 

In order to demonstrate the validity and correctness of the above derivations, we apply the 
new methods and the classical fourth-order Runge-Kutta method to a sample inhomogeneous 
linear scalar ODE given by 

— = Xu + ae^' (16) 

dt 


The exact solution is 


u(t ) = «( 0)e^‘ + 


p - A. 


(17) 


when p*A~ For the example, we use A, = — 1 +5 i, p = «, a - 10, and u(0)- 1. With these 
parameters, the exact solution at t = 2 is -2.60430984499756 - 0.20741 39 1939986*. In the table 
below, the magnitude of the errors obtained at t = 2 are compared for time step sizes ranging 
from 0.4 to 0.05. The classical fourth-order method is designated RK4, the new low-storage 
fourth-order method, RK4L, the new fifth-order method, RK5, and the new sixth-order method, 
RK6. In addition to showing the error magnitude, the table also shows the order of the method 
estimated from the error at the time step indicated and the next larger time step. The main 
results of the table are that the new methods approach the expected order of accuracy as the time 
step size is reduced and that the low-storage fourth-order method is as accurate as the classical 
method for inhomogeneous linear ODEs. 



RK4 

RK4L 

RK5 

RK 

6 

h 

error 

order 

error 

order 

error 

order 

error 

order 

0.4 

0.3437 

- 

0.3323 


0.2761 


0.0664 


0.2 

0.0343 

3.33 

0.0340 

3.29 1 

0.0059 

5.54 

8.63e-4 

6.27 

0.1 

0.0020 

4.11 

0.0020 

4.10 

1.70e-4 

5.12 

1.24e-5 

6.12 

0.05 

1.19e-4 

4.06 

1.19e-4 

4.05 

5.1 le-6 

5.06 

1.87e-7 

6.06 


Table 1. Sample computations using the new methods. 
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Conclusions 


Three new Runge-Kutta methods have been presented for the integration of linear systems of 
ODEs with constant coefficients. If the time step size is limited by stability, then the new 
fourth-order method is the most suitable of the new methods. This method requires less memory 
than the classical fourth-order Runge-Kutta method and less computational effort than the low- 
storage method proposed in [3]. If the time step is limited by accuracy, and memory is a 
secondary concern, then the new fifth- and sixth-order methods present an efficient new 
alternative. Since the expense of the methods is roughly proportional to the number of stages for 
the problems of interest here, the new fifth- and sixth-order methods are significantly more 
efficient than their counterparts for nonlinear ODEs. The sixth-order method is a particularly 
good choice for use with high-order spatial discretizations. 
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Figure 1: Stability contours for the fourth-order ( ), fifth-order ( ), and sixth- 

order (• • •) methods. 



coh/p 


Figure 2: Amplitude error produced by the fourth-order ( ), fifth-order ( ), and 

sixth-order (• ■ •) methods. 
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Figure 3: Phase error produced by the fourth-order ( ), fifth-order ( ), and 

sixth-order (■ • •) methods. 
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